Re-examining the electronic structure of germanium: A first-principle study 
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We report results from an efficient, robust, ab-initio method for self-consistent calculations of elec- 
tronic and structural properties of Ge. Our non-relativistic calculations employed a generalized 
gradient approximation (GGA) potential and the linear combination of atomic orbitals (LCAO) for- 
mahsm. The distinctive feature of our computations stem from the use of Bagayoko-Zhao-Williams- 
Ekuma-Frankhn (BZW-EF) method. Our results are in agreement with experimental ones where the 
latter are available. In particular, our theoretical, indirect band gap of 0.65 eV, at the experimental 
lattice constant of 5.66 A, is in excellent agreement with experiment. Our predicted, equilibrium 
lattice constant is 5.63 A, with a corresponding indirect band gap of 0.65 eV and a bulk modulus of 
80 GPa. We also calculated the effective masses in various directions with respect to the F point. 
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There has been great interest in germanium (Ge) due 
to the central role it plays in modern electronics. As 
a consequence of this interest, there has been many ex- 
perimental [T|-[T5] as well as theoretical [T21 fT^HTO] stud- 
ies, with the latter utilizing computational techniques 
of varying sophistications, ranging from tight binding 
method, the empirical pseudopotential method, and den- 
sity functional theory (DFT) methods to approaches that 
entirely go beyond DFT. 

Generally, DFT pUH^ is the most widely used ab 
initio theory in modern computational studies of elec- 
tronic properties of materials. Angle resolved photoe- 
mission and X-ray resonant emission spectra [U [5] show 
that the upper valence bands for traditional sp semicon- 
ductors and metals are well reproduced within the local 
density functional approximation (LDA) to DFT. How- 
ever, the low-laying conduction bands are severely under- 
estimated, leading to the well-known band gap problem. 

The Green function, G, and the screened Coulomb 
interaction, W, method (GW) has been used to accu- 
rately reproduce the band gap and other properties of 
Ge. However, the GW results do not lead to any sig- 
nificant changes in the dispersion of the upper valence 
bands jH] as compared to DFT findings, but rather a 
rigid shift in their energies and a reduction of the calcu- 
lated valence-band width [15]. In most GW methods, the 
Dyson equation is not solved self-consistently [16J , result- 
ing in a violation of charge conservation [53]. From the 
work of Ku and Eguiluz [TC] , it is evident (as pointed out 
by the authors) that even though their calculated band 
gap is in agreement with experiment, low-laying empty 
states near the X point are severely underestimated. 

Heyd-Scuseria-Ernzerhof (HSE) hybrid functional [571 
[55] calculations reported an indirect gap in the range 0.63 
- 0.77 eV [mSI]. The latter authors [^MH] did not pub- 
lish any electronic energies at the high symmetry points 
of Ge. The HSE approach involves a range separation 



of the exchange energy into some fractions of nonlocal 
Hartree-Fock exchange and of DFT exchange potentials. 
There are several versions of HSE ^M^M^^- The 
approach generally entails a range separation parameter, 
uj, which varies between and cxi. There exists a value 
of UJ that gives the correct gap for a given system; this 
value is adjusted from one system to another [301 EH [33] • 

The theoretical underestimations of band gaps and 
other energy eigenvalues have been ascribed to the in- 
adequacies of DFT potentials for the description of the 
ground state electronic properties of Ge [13 HH [IS]- 
Other methods [TH IH [THl [551 [301 135] that entirely go 
beyond DFT do not obtain the band gap value of Ge 
and related electronic structure quantities without ad- 
justments or fitting parameters. Further, to the best of 
our knowledge, there is no rigorous theoretical results 
predicting (or addressing) the indirect band gap of Ge. 
This unsatisfactory situation is a key motivation for our 
work. A concomitant motivation is to attain a highly ac- 
curate, ab-initio, self-consistent computational capability 
that lends itself to informing and guiding the design and 
fabrication of semiconductor based devices. 

The possibility for the referenced attainment has been 
stated by Bagayoko and co-workers [53H35] . despite per- 
ceived limitations of DFT in the literature. They ex- 
plained the fact that the derivative discontinuity of the 
exchange correlation energy has yet to be proven to be 
non-zero in semiconductors |36| . even though some be- 
lieve it to be non-zero. Perdew et al.^\, following a 
thought experiment on a diatomic molecule, established 
the existence of a derivative discontinuity of the exchange 
correlation energy. Generalizing this derivative disconti- 
nuity to the case of semiconductors, Perdew and Levy 
[30] showed that the exchange correlation potential may 
jump by the discontinuity, A^c, when the number of 
electrons in the system under study increases by one. 
In their work on this discontinuity in insulators. Sham 
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and Schliiter [JT] explicitly stated that while they es- 
tablished the existence of the discontinuity, they could 
not determine whether it is zero or not in real insula- 
tors. Subsequent work by Sham and Schliiter |42] de- 
rived the discontinuity of the functional derivative of Ei^c 
in insulators by considering an increase of the number 
of electrons by one. Cautiously, these authors concluded 
that the discrepancy between calculated and measured 
band gap is a measure of the discontinuity A^c - given 
the results from several calculations - if the employed 
LDA potentials are assumed to be good approximations. 
The description of our method below indicates the strong 
possibility that some current LDA and GGA potentials 
may be very good approximations. Computational ap- 
proaches may be sources of the discrepancy, as the per- 
ceived limitations of DFT, described in the literature, 
are far from being settled. This statement is supported 
in part by several DFT results, including predictions, in 
excellent agreement with experiment [36L 143) . An ampler 
discussion of the derivative discontinuity and of effects of 
computational approaches is provided in the supplemen- 
tary material. 

Ge has an experimental, indirect band gap (E^^'") of 
0.664 eV at room temperature [44]. Theoretical calcu- 
lations using several techniques have led to band gaps 
of Ge in the ranges of -0.02 to 0.35 eV for LDA and 
GGA [milS] and of 0.51 to 0.94 eV for the GW method 
[T51 [T5] . In most of these GW results, the position of 
the conduction band minimum {CBmin) is not at the L 
point. For instance, in the work of Ku and Eguiluz |16j . 
it is at the X point for their J^gwI'^gw], J^gwI'^lda], 
and J2gwI-^lda]+t^o 3d's methods since their predicted 
band gap E^~^ value of 0.71, 0.49, and 0.49 eV are 
smaller than the corresponding 0.79, 0.51, 0.51 eV for 
Eg^'". However, their PS-based + CPP method obtained 
correctly the position of CBj^in- A similar trend can be 
seen in the work of Fleszar [U] where E^"'" of 0.57 eV is 
smaller than E^^^ of 0.62 eV, and in cases where the or- 
dering is correct, the gap is over-estimated by ~ 0.3 eV. 
The work of Aulbur et al. [19 found the lowest eigen- 
value energies at the high symmetry points to be '-^ 1.0 
eV over-estimated. The same can be said of several other 
GW results [551 131] where the gap of Ge is found to be 
direct instead of indirect, as established by experiment. 

In this letter, we present a simple, robust approach 
based on basis set optimization [37j [38l [43] . We use this 
method to calculate the electronic and structural prop- 
erties of Ge and compare our results to experiments. 

In the ground state, Ge crystallizes in the diamond 
structure (space group: - Fd3m) [46ti48] with room 
temperature lattice constant of 5.66 A [471 . 

Our ab initio, self consistent, nonrelativistic calcula- 
tions employed a linear combination of atomic orbitals 
(LCAO). We utilize the electronic structure package from 
the Ames Laboratory of the US Department of Energy 
(DOE), Ames, Iowa [33] ■ For the LDA computations. 




FIG. 1: Calculated band structure of Ge, as obtained using 
our GGA BZW-EF optimal basis set. The horizontal, dot- 
ted line denotes the position of the Fermi energy (Ep) which 
has been set equal to zero. The calculated indirect band gap 
(E^""^) of 0.65 eV is basically the same as the room temper- 
ature experiment one of 0.66 eV [44] . 




FIG. 2: (Color online). Calculated density of states (DOS) 
and partial DOS of Ge, as obtained using our GGA BZW-EF 
optimal basis set. The vertical, straight, dashed line denotes 
the Fermi energy (Ep) which has been set equal to zero. As 
per the insert, it clearly shows the energy band gap. 

we used the Ceperley and Alder density functional ap- 
proximation (DFT) exchange-correlation functional fSH] 
as parameterized by Vosko-Wilk-Nusair [51]. We refer 
to it as the CA-VWN LDA potential. The GGA cal- 
culations were carried out using the Ceperley and Alder 
DFT exchange correlation contribution [50] as parame- 
terized by Perdew and Wang [551 - 151] . We refer to it as 
the CA-PW GGA potential. 

The distinctive feature of our calculations is the use 
of the Bagayoko, Zhao, Williams, Ekuma, and Franklin 
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(BZW-EF) method [351 |55]. In the BZW method, 
[571 1551 1151 [55H55] a minimum of three successive, self- 
consistent calculations are performed. For the first one, 
a relatively small basis set, no smaller than the mini- 
mum basis set that accounts for all the electrons in the 
system, is employed. This set is augmented with one or- 
bital for the second calculation. Depending on the s, p, 
d, or / nature of this orbital, 2, 6, 10, and 14 functions, 
respectively, are added to the basis set. The occupied 
energies of Calculation I and II are compared graphi- 
cally and numerically. They are generally different, with 
eigenvalues from Calculation II being lower or equal to 
corresponding ones from I. This process of adding one 
orbital and of performing self-consistent calculations is 
continued until a calculation is found, say N, to have the 
same occupied energies as Calculation (N-l-1) that im- 
mediately follows it. Then, the outputs of Calculation 
N provide the LDA/GGA BZW description of the ma- 
terial under study. The basis set for this calculation is 
the optimal basis set, i.e., the smallest one that leads to 
the minima of all the occupied energies. In Calculation 
(N-M) and others with larger basis sets that include the 
optimal basis set, the charge density, the potential, and 
the Hamiltonian do not change as compared to their val- 
ues from Calculation N; nor do the occupied energies. 
These calculations, however, do lead to the lowering of 
some unoccupied energies on account of the basis set and 
variational effect stemming from the Rayleigh theorem 
[551 jSD] . For Ge, the basis set of Calculation III was the 
optimal one [Ge {3s3p3d4:s4p4:d5p) , for the valence states 
only], where the Ad and 5p orbitals are unoccupied. 

The enhancement of the original BZW method is in the 
methodical increase of the basis set in our calculations. 
This enhancement leads to adding p, d, and / polar- 
ization orbitals, for a given principal quantum number, 
before adding the spherically symmetric s orbital. These 
additional unoccupied orbitals are needed to accommo- 
date the reorganization of the electron cloud, including 
possible polarization, in the crystal environment. For 
valence electrons in molecules to solids, polarization has 
primacy over spherical symmetry. The BZW method has 
been shown to lead to accurate ground state properties 
of many semiconductors: c-InN [59j, w-InN ^58^, w-CdS 
[57], c-CdS [SB, rutile-Ti02 [13], SrTi03[55] AlAs [52], 
GaN, Si, C, Ru02 [S5^, and carbon nanotubes [35]. De- 
tails of our method have been explicitly explained in the 
Uteratures [33 [Ml 113 113 |56l IfiQl 115 ■ 

The Brillouin zone (BZ) integration for the charge den- 
sity in the self consistent procedure is based on 28 special 
k points in the irreducible Brillouin zone (IBZ). The com- 
putational error for the valence charge is 2.3 x 10~^ eV 
per valence electron. The self consistent potential con- 
verged to a difference of 10~^ after several tenths of it- 
erations. The energy eigenvalues and eigenfunctions are 
then obtained at 161 special k points in the IBZ for the 
band structure. A total of 89 weighted k points, cho- 
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FIG. 3: (Color online). 
Calculated total energy 
per unit cell (Et) as a 
function of the lattice 
constant (A) of Ge, as 
obtained with the GGA 
BZW-EF optimal basis 
set. The calculated 
equilibrium lattice con- 
stant is 5.63 A, is little 
smaller than the room 
experimental value of 
5.66 A, as expected. 



sen along the high symmetry lines in the IBZ of Ge, are 
used to solve for the energy eigenvalues from which the 
electron density of states (DOS) are calculated using the 
linear analytical tetrahedron method [BB]. The partial 
density of states (pDOS) and the effective charge at each 
atomic sites are evaluated using the MuUiken charge anal- 
ysis procedure [67]. We also calculated the equilibrium 
lattice constant Uo, the bulk modulus (Bo), the associated 
total energy, and the electron and hole effective masses 
in different directions. 

In calculating the lattice constant, we utilized a least 
square fit of our data to the Murnaghan's equation of 
state [S5]. The lattice constant for the minimum total 
energy is the equilibrium one. The bulk modulus (Bo) is 
calculated at the equilibrium lattice constant. 

The difference in energy between the GGA BZW-EF 
and LDA BZW-EF computational results is small ^ 0.05 
eV. Thus, the results reported here are those from GGA 
BZW-EF unless otherwise stated. The results of the elec- 
tronic structure are given in Figs. (Tj and [2| Total energy 
versus lattice constant data are as shown in Fig. [3] while 
Table[l|shows the eigenvalue energies at various high sym- 
metry points in the Brillouin zone. 

Our ab-initio, first principle results show that the fun- 
damental gap of Ge is an indirect one with the maximum 
of the valence band (VBmax) occurring at F and the min- 
imum of the conduction band {CBmin) at the L point (cf. 
Figs, [ijand [2] and Table |l]). Our calculated indirect gap, 
using the equilibrium lattice constant, is 0.649 eV, while 
the indirect band gap value using the experimental lat- 
tice constant is 0.652 eV. Our calculated direct gap at 
the r point is 0.857 eV (cf. Table |l|. Our calculated in- 
direct band gap value using the LDA BZW-EF is 0.644 
eV with a corresponding direct gap of 0.857 eV. The LDA 
and GGA potentials lead to the same direct gap, while 
the LDA indirect gap is 8 meV smaller than that for 
GGA. Our calculated values are in good agreement with 
experimental room temperature indirect band gap value 
of 0.664 eV and direct band gap value in the range 0.805 - 
0.895 eV [111169]. The accurate, calculated values for the 
indirect and direct gaps of Ge are a significant improve- 
ment over some previous results that employed methods 
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TABLE I: The calculated electronic energies (eV) of different bands at high symmetry points in the Brillouin zone of Ge, as 
obtained by GGA BZW-EF calculations, compared to experimental ones where the latter are available. Experimental data are 
from Refs. [31 [H iZl [Ml ES] • 
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(or approaches) beyond DFT [Iil[Tni[ini[251[3S]. Our cal- 
culated spectra shown in Fig. [2] especially for the valence 
states, are in good agreement with experiments [31 (H |6l- 
[Hl [70]. The same can be said of our calculated bands, as 
shown in Fig. [T] which are in good agreement with exper- 
iment [51 [7T] . The occupied bandwith is calculated to be 
12.56 eV, in very good agreement with room temperature 
values of 12.6 ± 0.3 |3l [3 |47l [64]. 

From Fig. [2j in the DOS of valence states, we predict 
two weak shoulders at -0.31 and -0.67 eV, and a peak at - 
1.50 eV. This peak is followed by relatively broad ones at 
-1.90, -2.41, and -2.86 eV, before two sharp peaks at -4.06 
and -7.41 eV. In the conduction bands, a small shoulder 
can be found at 0.71 eV followed by a peak at 2.10 eV. 
Other peaks can be seen at 3.07, 3.88, 4.44, and 5.06 eV. 
The position of these critical points are in good agree- 
ment with experiments [SHSl [El [Ml IMl llSl Il3] . For in- 
stance, scanning tunneling spectroscopy measurements of 
Feenstra [5] reported peaks in the valence states at -0.47 
± 0.06 eV and -1.23 ± 0.030 eV while the inverse pho- 
toemission spectra of Jackson et al. [73] reported peaks 
at 2.40 ± 0.10, 4.20 ± 0.10, and 5.50* ± 0.10 eV. We 
note that our spectra are not broadened as it is the case 
in (many) experiments; they may therefore have more 
features. 

From Fig. [2] significant contributions of 3d states can 
be seen in both the valence and conduction bands. Even 
though Ge is known to be an sp material, the correct 
treatment of 3d states and others as valence ones is crit- 
ical for accurate results. Crystal symmetry gives rise 
to polarizations that are absent in isolated atom. In our 
work, only Ge (Is2s2p) were considered to be in the core. 



The importance of this view has been highlighted by the 
work of Shirley et al. [TJ] and confirmed by those of Ku 
and Eguiluz [16] and of Rohlfing et al. [25] . It generally 
leads to a better description of the band structure en- 
ergies, especially in obtaining the fundamental, indirect 
gap of Ge. 

Effective masses are measures of the accuracy of a 
method as they depend very sensitively on the curva- 
tures of the bands. They determine in part the transport 
properties, the Seebeck coefficient, and the electrical con- 
ductivity of materials. The calculated electron effective 
masses at the bottom of the conduction band along the F 

- L,T - X, and T - K directions, are 0.043 - 0.051, 0.043 

- 0.051, and 0.043 - 0.052 (all in units of the electron 
mass), respectively. These results are in good agreement 
with experiment |47j . 

The total energy versus the lattice constant data are 
shown in Fig. [3] The data fit well to the Murnaghan 
equation of state (EOS) j68j. The calculated bulk modu- 
lus. Bo, is 80 GPa, in close agreement with experimental 
values of 75 - 80 GPa [HIZSlTS]. We also calculated the 
equilibrium lattice constant from the total energy min- 
imization. The calculated equilibrium lattice constant 
is 5.63 A. The experimentally reported room tempera- 
ture lattice constants are in the range of 5.62 to 5.66 A 

dUiS]. 

In summary, our results show that the electronic and 
related properties of Ge can be accurately described, with 
LDA/GGA potentials, by a careful search for a basis set 
that is verified to be complete for the description of the 
ground states. 
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SUPPLEMENTARY MATERIAL 

Despite the great progress made possible by density functional theory (DFT), from 1964 to present, problems 
associated with obtaining theoretically the measured energy or band gaps, for finite and crystalline semiconductors, 
respectively, have persisted. Specifically, most DFT calculations, with emphasis on those utilizing local density 
approximation (LDA) and semi-local potentials, have led to semiconductor band gaps that are 30 - 50% smaller than 
their corresponding, measured values. Much effort has been deployed to find explanations of and remedies to this 
recalcitrant band gap problem. Perdew and Zunger [77j introduced the self interaction correction (SIC) to local spin 
density (LSD) approximation calculations. While the exact functional for the ground state is self interaction free, 
these authors discussed corrections that appear to be needed for the description of finite systems, beginning with 
atoms, and of localized states in solids. This self interaction is argued to contribute to the underestimation of the 
band gaps of insulators by DFT calculations jTTj. Consequently, self interaction corrections (SIC) are expected to 
improve the agreement between calculated band gaps and measured ones, in addition to improving binding energies 
and bringing orbital energies closer to removal energies [77J [TS] . While self interaction corrections have led to some 
improvements in band gap calculations, they have not totally resolved the problem. Applications of SIC have mostly 
overestimated the band gap of semiconductors ^79l [80] . According to Cohen et al. [8T] , self interaction is well-defined 
only for one-electron systems. 

According to the literature, a major source of the theoretical underestimation of band gaps consists of the derivative 
discontinuity of the exchange correlation energy, Exc |39fl42j . Perdew et aZ.[39j, following a thought experiment on 
a diatomic molecule, established the existence of a derivative discontinuity of the exchange correlation energy, i.e., a 
discontinuity in the exchange correlation potential, Yxc- Perdew and Levy [40! generalized this discontinuity to the 
case of semiconductors. They showed that the exchange correlation potential may jump by the discontinuity, A^;^, 
when the number of electrons in the system under study increases by one. Band gaps calculated with a local density 
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approximation (LDA) potentials, according to their findings, are to be augmented by this discontinuity in order to 
reproduce the corresponding, measured values. The authors suggested, without claiming to have a proof of it, that 
this discontinuity is a non zero (and positive) in real semiconductors and insulators. Sham and Schliiter |41) also 
found a derivative discontinuity of E^c in semiconductors. These authors, however, asserted that their work does not 
show whether or not this discontinuity is non zero in real insulators. Subsequent work by Sham and Schliiter |42| 
derived the discontinuity of the functional derivative of E^c in insulators by considering an increase of the number of 
electrons by one. Cautiously, these authors concluded that the discrepancy between calculated and measured band 
gap is a measure of the discontinuity A^c - given the results from several calculations - if the employed LDA potentials 
are assumed to be good approximations. The description of our method below indicates the strong possibility that 
some current LDA and GGA potentials may be very good approximations. 

Despite its popular use to explain the disagreement between calculated and measured band gaps, the above discon- 
tinuity has not yet been established to be non zero in real semiconductors or insulators. Further, Sham and Schliiter 
[42j underscored the fact that, in principle, DFT and Kohn Sham LDA hold only if the number of particle is kept 
constant. The question could arise whether or not the discontinuity, derived by considering a change of the number of 
particle, is strictly applicable to DFT or LDA calculations. From the preceding, it has not yet been established that 
DFT or LDA calculations cannot obtain the correct band gaps, despite the fact that presently known LDA potentials 
do not have a discontinuity and that most of the numerous, previous ab-initio DFT and LDA calculations did not. 

Another presumed contributor to the band gap underestimation by theory stems from the use of local (LDA) and 
semi-local (GGA) potentials. The question naturally arises as to what extent the local and semi-local potentials fail 
to capture key feature of the exact one. We are aware of no definitive answer, given that the exact one is not known. 
We would have had to delve into this matter further if we were dealing with molecules or their dissociation. The solid 
state systems of interest to us, to judge by previous results obtained with our method [351 [SS], possible errors due to 
the use of local and semi-local potentials appear to be very small. 

There exist several approaches that have been introduced to address the band gap problem. Review articles and 
books are the best sources for discussing these approaches and for examples of the many DFT calculations that led to 
band gaps much smaller than their corresponding, experimental counterparts. In contrast, a summary of results from 
BZW LDA calculations for over 10 materials show agreement between theory and experiment. Illustrative examples of 
discrepancies between theory and experiment follow. The case of Ge is summarized in this article. Some previous LDA, 
GGA, and GW calculations did not yield the measured band gap, from first principle. A table provided by Ekuma 
and Bagayoko P3| shows a multitude of DFT calculations with vastly different band gaps for titanium dioxide. With 
the computational method described here, Ekuma and Bagayoko obtained the measured, direct gap and predicted an 
indirect one. For elemental silicon, Griining et aL[82] reported an LDA band gap of 0.7 eV, much smaller than the 1.25 
eV they reported as the measured value. These authors also performed calculations with the exact exchange (EXX), 
EXX plus LDA, EXX plus the random phase approximation (RPA). The last approach or scheme yielded 0.6 eV, a 
gap smaller the above LDA gap, while the first two led to 1.5 and 1.6 eV respectively, values much larger than the 
experimental one. With the original version of our method, Zhao et al. utilized an LDA potential to obtain a gap of 
1.02 for Si, much closer to the experimental one. Generalized gradient approximation (GGA) calculations have led to 
improvements of calculated properties of materials, including lattice parameters. Specifically, Hao et a/. [83] reported 
revised Tao-Perdew-Staroverov-Scuseria (revTPSS) meta-GGA calculated lattice parameters that are in agreement 
with experiment following a zero-point phonon correction, for over 50 materials. Despite this very significant success, 
most GGA and meta-GGA calculations, including the previous ones discussed here for Ge, have not produced band 
gaps in agreement with experiment. 

From the above summary, historical overview of the band gap problem, it appears that the scientific community 
believes that the derivative discontinuity of the exchange correlation energy is the main source of the disagreement 
between DFT calculated energy and band gaps and their corresponding, measured ones. This belief led to the 
development of several schemes aimed at resolving the band gap problem. Except for the few, most of these schemes 
are ad hoc as they include adjustable parameters that vary with the material under study. The continuing growth in 
the number of these schemes seems be a problem in itself, the ad ho nature of most of them does not lend itself to 
predictive capabilities from first principle, the aim of theory to inform and to guide experiment. The only exception 
of the above trend consists of the work of our group. This work has not yet gotten the attention of the community 
at large, presumably due to the strength of the above belief, one the one hand, and the preponderance of results 
that are explained with the discontinuity. As we previously noted [58j, the situation resembles that of the Ptolemaic 
model of the solar system where epicycles were continually introduced to explain its disagreement with observations. 
The quintessential point in support of the our method, described below, is the following: For all DFT calculations 
of energy bands, the minima of the occupied energies, which add up to yield the ground state energy of the electron 
system, are obtained from the theory if the "correct" ground state charge density is utilized, subject to the constraint 
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that the number of particle is kept constant |201 121] . 

Most of the previous DFT calculations, including those with GGA and LDA potentials, have consisted of judicious 
selecting large basis set and of performing iterations to obtain self consistent eigenvalues of the Kohn-Sham type 
equation. It is assumed that the single basis set in question leads to the correct representation of the electronic cloud 
in the system under study, a system that can be drastically different from an atomic or ionic one. In particular, as we 
recently pointed out, polarization (p, d, and f orhitals) has primacy over spherical symmetry (s orbital) for systems 
varying from diatomic molecules to solids. Hence, utilizing basis sets derived from calculations of properties of atoms 
for the study of solids is potentially problematic. Indeed, the angular symmetries in these systems are vastly different 
from those for atoms and ions. The need for the method described below becomes apparent with the realization 
that, irrespective of the degree of convergence of the iterations, a single trial basis set that has a major symmetry 
inadequacy for the description is not going to lead to physically valid DFT eigenvalues as the implacable condition of 
using the "correct" ground state density will not be met. 

The original version of our method, named after Bagayoko, Zhao, and Williams (BZW) was introduced in 1998 [SD] 
and further explained in 1999 [SHj. The method consists of implementing the linear combination of atomic orbitals 
(LCAO) formalism by methodically searching for the smallest basis set, called the optimal basis set, that leads to the 
minima of the occupied energies. This search begins with a deliberately small basis set that is not smaller than the 
minimum basis set, i.e., the smallest one needed to account for all the electrons in the system. The self consistent 
calculation with this basis set is followed by another whose basis set uses the previous one plus one additional orbital. 
The occupied energies from the two calculations are compared numerically and graphically. In the more than 20 
systems we have studied, these occupied energies from these first two calculations have been different, with those of 
calculation II generally lower than the corresponding one from Calculation I. Calculation III is then carried out, using 
the basis set in Calculation II plus an additional orbital. The occupied, self consistent energies from Calculations II 
and III are also compared. This process of augmenting a basis set, performing a new calculations, and comparing its 
results with those of the one immediately preceding it, continues until a calculation is found, say N, to have exactly the 
same occupied energies as the one immediately following it. This perfect superposition establishes that the minima 
of the occupied energies have been reached and that the corresponding basis set give the best representation of the 
ground state charge density of the system. 

Before elaborating further on the physical content of the method, we note that adding an orbital means increasing 
the size of the basis set (and hence the dimension of the Hamiltonian) by 2, 6, 10, or 14 depending on the s, p, 
d, or / character, respectively, of the orbital in question. In the original BZW, we added orbitals in the order of 
their energies resulting from the atomic calculations, i.e., the orbitals corresponding to the lowest laying excited, 
atomic state were successively added. As apparent from our earlier results, the BZW method practically led to the 
measured band gap of semiconductors we studied. Further, our predictions of the band gaps and other properties for 
cubic Si3N4 [53] and cubic InN jSS] were totally confirmed by experiment [55H57] . Following the works of Ekuma and 
Franklin [371 [Ml 1331 ISI] i we realized that valence electrons in multi-atomic systems simply do not follow the symmetry 
landscape that prevail for isolated atoms or ions. The aim is to obtain a better representation of the electronic cloud 
(ground state charge density) of the system under study. Hence, in most of our subsequent calculations, for a given 
principal quantum number, we add p, d, and / (if applicable) orbitals before the s orbital for that principal quantum 
number. This counter-intuitive ordering, for isolated atoms, is simply needed for 'some' multi-atomic systems. The 
initials of Ekuma and Franklin (EF) are added to the name of the enhanced method [BZW-EF] in recognition of their 
extensive calculations whose results led Bagayoko to see the necessity for this new order. While the BZW method 
led to band gaps that were sometimes smaller by 0.1 - 0.3 eV, insignificant differences between BZW-EF calculated 
gaps and corresponding experiment ones are in the second decimal place for the systems studied to date. The BZW- 
EF method led to an upper valence band width more than 1.0 eV larger than was obtained with the BZW, with a 
significant improvement in agreement with experiment. 

The origin of the changes in the band structure and the band gap when the basis set increases toward the optimal 
one consists of the progressively better representation of the ground state charge density. As per the derivation of 
DFT, the minima of the occupied energies are obtained if "correct" charge density for the ground state is employed. 
These changes are due to physical interactions, given that the charge density, the potential and hence the Hamiltonian 
change from one calculation to the next. We should underscore here that while our focus is on occupied energies (i.e., 
DFT is a fundamentally ground state theory) , when these energies reach their minima, so do the low laying unoccupied 
energies, up to 9 - 10 eV for the materials studied to date with the BZW-EF method. For the many systems with the 
BZW, most low laying unoccupied energies also converged up to 5 - 6 eV, as was the case for wurtzite indium nitride 
[88j . For a few materials, this convergence of the lowest unoccupied energies was not achieved with the original BZW 
method. As shown by Zhao et aL[56_ the low-laying unoccupied energies converge when the occupied ones do, due to 
the fact that at least one band crosses the Fermi level. This fact may partly explain the early successes of DFT in 
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describing metals as compared to semiconductors. 

The description of what occurs when basis sets much larger than the optimal ones are employed will complete the 
description of the BZW-EF method. We first recall that the basis set immediately following the optimal one leads 
to the occupied energies obtained with the optimal one and to the same unoccupied energies up to 9 - 10 eV. So, 
by much larger basis sets, we mean the ones that are larger than the basis sets immediately following the optimal 
ones. Earlier works by Zhao and Co-workers |561 160| verified that basis sets larger than the optimal one do not change 
the charge density, the potential, and the Hamiltonian, nor do they change the occupied energies. In the absence of 
changes in the Hamiltonian, i.e., the physics of the study, the additional lowering of unoccupied energies with these 
much larger basis set cannot be ascribed to DFT. However, the Rayleigh theorem provides an explanation of the 
unphysical lowering of these energies. The theorem states that when the same eigenvalue equation is solved with two 
basis sets of different size, such that the larger one includes the smaller one, then the eigenvalue obtained with the 
larger basis set are lower or equal to their corresponding one obtained with the smaller basis set. 

Clearly, after the optimal basis set is reached, and that the occupied energies are no longer changed from their 
values obtained with the optimal one, the lowering of unoccupied energies can be ascribed to a mathematical artifact 
that is the manifestation of the above theorem. We therefore contend that this extra lowering, a variational basis 
set effect [SHI HOj: is a major source of discrepancies between many previous DFT calculations and between these 
calculations and experiment, as far as band gaps are concerned. The lowest laying conduction bands, with full physical 
meaning when they result from the optimal basis set, partly lose their physical content due to the above effect. It is 
important to note that this is the case for any basis set that is not the optimal one, whether it is smaller or larger 
than the optimal, or simply lacks orbital or angular features of the optimal one. For basis sets that do not totally 
include their corresponding optimal ones, even the occupied energies are not totally DFT results. The preceding 
lines in this paragraph point to the great difficulty in obtaining physically meaning DFT occupied and low laying 
unoccupied energies with a single trial basis set, irrespective of the degree of convergence of the applicable iterations 
for self consistency. 



